Exact solution for quantum dynamics of a periodically-driven two-level-system 
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We present a family of exact analytic solutions for non-linear quantum dynamics of a two-level 
system (TLS) subject to a periodic-in-time external field. In constructing the exactly solvable 
models, we use a "reverse engineering" approach where the form of external perturbation is chosen 
to preserve an integrability constraint, which yields a single non- linear differential equation for the 
ac- field. A solution to this equation is expressed in terms of Jacobi elliptic functions with three 
independent parameters that allows one to choose the frequency, average value, and amplitude of 
the time-dependent field at will. This form of the ac-drive is especially relevant to the problem of 
dynamics of TLS charge defects that cause dielectric losses in superconducting qubits. We apply 
our exact results to analyze non-linear dielectric response of such TLSs and show that the position 
of the resonance peak in the spectrum of the relevant correlation function is determined by the 
quantum-mechanical phase accumulated by the TLS wave-function over a time evolution cycle. It 
is shown that in the non-linear regime, this resonance frequency may be shifted strongly from the 
value predicted by the canonical TLS model. We also analyze the "spin" survival probability in the 
regime of strong external drive and recover a coherent destruction of tunneling phenomenon within 
our family of exact solutions, which manifests itself as a strong suppression of "spin- flip" processes 
and suggests that such non-linear dynamics in LC-resonators may lead to lower losses. 

PACS numbers: 03.65.Yz, 42.50.Hz, 03.67.Lx 



I. INTRODUCTION 

The problem of a periodically-driven two-level system 
(TLS) appears in many physical contexts including mag- 
netism, superconductivity, structural glasses and quan- 
tum information theoryi - — The interest in this old prob- 
lem has been revived recently due to advances in the field 
of quantum computing (see, e.g., Refs. j8Ml2J and refer- 
ences therein). First of all, a qubit itself is a two- level 
system and the question of its evolution under an exter- 
nal time-dependent perturbation is obviously of interest. 
Also, the physical mechanism that currently limits coher- 
ence particularly in superconducting qubits is believed to 
be due to other types of unwanted TLSs within the qubit, 
whose charge dynamics under a periodic-in-time electric 
field gives rise to dielectric losses directly probed in exper- 
iment^^ In what follows, we mostly apply our solution 
to the latter charge TLS model, but the general methods 
and some particular results of this work evidently can be 
applied to a much broader range of problems (see, e.g., 
Ref. [15j and references therein). 

One of the key metrics of a superconducting qubit is 
the quality factor, which is defined as a ratio of the real 
and imaginary parts of the dielectric response function, 
e(oj), evaluated at the resonant frequency of the corre- 
sponding LC-circuit, Q = Kes(uj T )/lme(uj r ). Very high 
values of the quality factor are required for the qubit to 
be operational. However, existing experiments consis- 
tently show significant dielectric losses that occur in an 
amorphous dielectric (e.g., in AI2 O3) used as a barrier in 
the Josephson junctions. It is believed that the losses are 
primarily due to the presence of charge two-level system 
defects in the barrier and/or the contact interfaces, which 
respond to an AC electric field in the LC-resonator. It 



is still unclear what the physical origin of these defects 
is, but an early work of Phillips 16 as well as very recent 
comprehensive density functional theory studies of Mus- 
gravei£ point to the OH-rotor defects as a very likely 
source of the dielectric loss. To determine the physical 
origin and the properties of the TLSs responsible for the 
dielectric loss is one of the central questions in the field 
of superconducting quantum computing and it has been 
larger ly the main physical motivation for our work. 

The usual theoretical approach to calculating the qual- 
ity factor and more generally the full dielectric response 
function, e{uo), involves a formal mapping of charge dy- 
namics in a double-well potential onto the problem of 
"spin" dynamics in an AC field, described by the "spin" 
Hamiltonian H(t) = h(t)-&/2, where & denotes the Pauli 

matrices and b(t) = 2 ( A^,0,e + <ixLS • E(t) ) is an ef- 
fective "magnetic field" that drives TLSs, with e, A t , 
and (ixLS being the TLS energy splitting, the tunneling 
amplitude between its two states, and the TLS dielec- 
tric moment correspondingly and E(t) is the AC electric 
field. A linear analysis within the canonical TLS model 
predicts that the dielectric function due to identical TLSs 
is peaked at the frequency, v = \/ A 2 + s 2 . Ad-hoc inclu- 
sion of T\ and T2 relaxation processes and the assump- 
tion about random distribution of TLS energy-splittings 
and tunnelings (typically assumed to be uniform and 
log-uniform correspondingly) lead to the quality factor 
Q oc yl + (Eo/E c ) x , with x ~ 2, Eq being the ampli- 
tude of an applied AC electric field and E c is a critical 
value of the amplitude which also encodes the informa- 
tion on the strength of the relaxation processes (see, e.g., 
Ref. [5]). Both formulas are used widely in interpreting 
experimental data and probing energetics of the relevant 
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FIG. 1: Schematic representation of an OH-rotor two- level 
system in an AI2O3 oxide . 16 i 17 Here, the role of the general- 
ized variable is assigned to the angle defined as an angle 
between the OH-bond and an axis perpendicular to the verti- 
cal A10 bond. At low enough temperatures, the phase space 
an isolated rotor is reduced to the two-states corresponding to 
the minima of the double- well potential V(9). Application of 
external ac- field parametrically coupled to the rotor's dipole 
moment induces oscillations between the two minima. 



TLS defects^i While this linear analysis is a fine ap- 
proximation to describe a majority of regimes currently 
studied experimentally, the existing experiments are cer- 
tainly capable and some dc3 access non-linear regimes as 
well, where the energy of the applied electric field is com- 
parable or larger than the relevant TLS energies. Hence, 
this non-perturbative regime is of clear experimental and 
theoretical interest. More importantly studies of non- 
linear dynamics may provide another effective means to 
probe the properties of TLSs. 

The mathematical formulation of the non-linear TLS 
dynamics problem studied in this paper is deceptively 
simple: We wish to solve the Schrodinger equation for 
a spinor wave-function, ^ = (^ + ), id t ^ = |b(t) • cr \I/ 
that describes a half-integer spin subject to a periodic- 
in-time magnetic field of the form, b(t) = 2 (At, 0, /(£)), 
where A t is a constant describing the coupling between 
the two states and the function f(t) = f(t + T) describes 
the time dependent perturbation. Despite the simplic- 
ity of the formulation, the problem is generally unsolv- 
able in analytic form for most cases of practical interest. 
The origin of this surprising fact can be understood if 
we introduce a new function R(t) = ip+(t)/ijj-(t), which 
reduces the matrix Schrodinger equation to the Riccatti 
equation d(_i t \R = 2fR + A t [l — R 2 ] . It is a non-linear 
differential equation that has known analytic solution in 
a very limited number of cases (note that the case of a 
monochromatic perturbation is not one of them). There- 
fore, to solve for TLS dynamics driven by a specific non- 
equilibrium field is equivalent to generating a particular 
solution to the Riccatti equation corresponding to this 



perturbation. Clearly this is a challenging mathematical 
task and this observation partially explains the current 
deficit of exact mathematical results. The difficulties 
in obtaining exact solutions have led to the emergence 
of several perturbative approaches, used in particular to 
characterize relaxation and dephasing rates in qubits as 
a function of driving amplitude (see, e.g., Ref. [15| and 
references therein). These analyses provide very useful 
physical insights and correctly describe the physics if the 
time-dependent perturbation is weak, but it is also clear 
that there exist non-linear effects beyond perturbation 
theory and it is desirable to have exact results to access 
this qualitatively different physics. 

The mathematical approach that we use in this paper 
to obtain exact results is to "reverse engineer" exactly 
solvable Hamiltonians of specific form relevant to the 
problem of interest. A key observation in our analysis is 
that finding a Hamiltonian corresponding to a given solu- 
tion is much easier than solving the Schrodinger equation 
with a given Hamiltonian. In some generalized sense, the 
two procedures are related to one another much like dif- 
ferentiation relates to integration. To see this, it is use- 
ful to consider the evolution operator, or the ^-matrix, 
which relates the initial state at t — to a final state 
at t > as follows, #(£) = S(t)V(0). In the absence 
of relaxation process the time-evolution is unitary and it 
satisfies the Schrodinger equation, idtS(t) = H(t)S(t). If 
we choose an arbitrary S-matrix, S = exp [— 1<&(£) • <r] G 
5/7(2)2, we can immediately reconstruct the correspond- 
ing Hamiltonian that gives rise to such evolution as fol- 
lows %{t) = id t S(t)S^(t). Using this method, one can 
generate an infinite number of exact non-equilibrium so- 
lutions and explicit models. These solutions may be of 
importance to physics of NMR, to the question of physi- 
cal implementation of gate operations on a qubit as well 
as of some mathematical interest. Nevertheless without 
additional constraints such analyses would generally pro- 
duce Hamiltonians of little importance to the problem of 
dynamics of TLS charge defects. 

A very useful insight that allows us to construc- 
tively narrow down the range of relevant dynamical sys- 
tems comes from the mathematically related problem of 
far- from- equilibrium superconductivity^"—. It is well- 
known that the reduced BCS Hamiltonian is algebraically 
equivalent to an interacting XY-spin model in an effective 
"inhomogeneous" magnetic field in the z-direction, whose 
profile is dictated by the bare single particle-energy dis- 
persion. Far from equilibrium, dynamics of a given An- 
derson pseudospin 22 is determined by an effective time- 
dependent self-consistent field of other pseudo-spins that 
it interacts with^ In many cases (determined by specific 
initial conditions), this BCS self-consistency constraint 
dynamically selects a specific order-parameter, such that 
the dynamics of essentially infinite number of spins is 
equivalent to the dynamics of few spins only.— For special 
sets of initial conditions, these spins move in unison and 
therefore the self-consistent "magnetic field" (or super- 
conducting order parameter in the language of BCS the- 



ory) is periodic in time. The reduced BCS model is inte- 
grable and there exists a very elegant prescription for con- 
structing exact non-equilibrium solutions to it, developed 
primarily by Yuzbashyan and collaborators i 2Q i 21 i 24 These 
solutions contain, in particular, exact spin dynamics in 
a periodic time-dependent field that can be expressed in 
terms of elliptic functions. In this paper, we generalize 
such anomalous soliton solutions of Yuzbashyan 25 to en- 
compass a wider range of time dependencies relevant to 
the problem of TLS dynamics, which is of our primary 
interest. 

This paper is organized as follows: Sec. II summa- 
rizes a general mathematical structure behind the "re- 
verse engineering" approach to constructing exact solu- 
tions for non-linear TLS dynamics. The specific Ansatz 
and technical details of our particular family of solu- 
tions for periodically-driven TLS dynamics are given in 
Sec. III. In Sec. IV, we use some representative solutions 
to illustrate the emergence of the coherent destruction 
of tunneling phenomenon. We also derive the spectrum 
of exact dielectric response function due to an ensemble 
of identical charge TLS in the presence of dissipation, 
which is introduced phenomenologically. In Sec. V we 
provide a summary of our results. In the Appendices we 
list some technical details of our calculations as well as 
useful relations aimed to shed more light on the subtle 
features of our theory. 



II. GENERAL FRAMEWORK FOR 
CONSTRUCTING EXACT SOLUTIONS 

In this paper, we derive a family of exact solutions for 
the non-dissipative TLS dynamics subject to an external 
ac-field. The main ingredient of our approach is a spe- 
cial Ansatz for the TLS's dynamics that corresponds to 
periodic-in-time but non-monochromatic external fields. 
Before proceeding to the specific Ansatz, let us first in- 
troduce a general algebraic framework for "reverse engi- 
neering" of exact solutions. We are interested in solving 
the non-equilibirum Schrodinger equation for the spinor 



idt*(t) = u(t)*(t), 



•w = (£ 



(i) 



where the Hamiltonian isH(i) = |b(£)-<r. As mentioned 
in the introduction, instead of solving Eq. (pQ) for the 
wave-function, we can consider the Schrodinger equation 
for the evolution operator that relates the initial and final 
states, #(£) = S(t)^(0). This equation for the ^-matrix 
has the form identical to Eq. (pp): 



id t S(t) = H(t)S(t), and S(0) = 1 



(2) 



but now it is an equation for the matrix function 
*§(£), which belongs to the two-dimensional represen- 
tation of the SU(2) group, while the Hamiltonian ex- 
pressed in terms of SU '(2)2 generators belongs to the two- 
dimensional representation of the su(2) algebra. Note 



that the form of Eq. (J2j) is such that it may be general- 
ized to an arbitrary spin or equivalently to an arbitrary- 
dimensional representation of SU{2) or it can be viewed 
as an equation of motion in the abstract group such 
that 7iabs0) = b W • ^abs e su(2) and S ahs (t) = 
exp —i<fr(t) -Jabs € 5/7(2), where J a b s are the corre- 
sponding generators. Therefore, a solution of the prob- 
lem in a particular representation, i.e., an explicit form 
of &(i), immediately gives the corresponding solutions 
in all other representations (e.g., a two- level- system dy- 
namics uniquely determines a "<i-level system" dynamics 
in the same field). This TLS problem that we are inter- 
ested in corresponds to the two-dimensional generators 

J a — \®ol with a a (a = x,y,z) being the familiar 
Pauli matrices. 

The problem of determining the solution, <£(£), from 
the magnetic field time-dependence h(t) is a complicated 
one, but the inverse problem is almost trivial. Indeed, 
if we select a specific S'-matrix (defined uniquely by the 
choice of a specific function, $(£)), the Hamiltonian will 
read 



H(t) = idtS(t)&(t), 



where 



S(t) = exp 



--*(*)•* 



(3) 



(4) 



Using the algebraic identities for the Pauli matrices, we 
obtain the corresponding magnetic field 

b(t) = i> n + sin $ h + (1 - cos $) [n x h] , (5) 

where *(£) = $(t)n(t), with \n(t)\ = 1. Note that 
one can generate exactly-solvable models by simply pick- 
ing an arbitrary &(t) dependence and using Eq. (J3j) to 
find the corresponding Hamiltonian. However, without 
guidance or luck, such an analysis would generally pro- 
duce complicated non-equilibrium fields that have lit- 
tle to do with an underlying physical problem. Let us 
however mention here that this procedure may be of 
interest to quantum computing in general, because the 
time-evolution governed by an ^-matrix can be viewed 
as a "gate operation" on the spin (if the TLS/spin cor- 
responds to a qubit rather than to a defect within a 
qubit). By picking "trajectories," 3>(£), on the algebra 
that start in the origin, i.e. 3?(0) = 0, but end at a spe- 
cific point at a time T, one can immediately determine 
the non-equilibrium magnetic pulse, b(t), or a class of 
such pulses, that will give rise to a desired gate operator 
G = S(T) = exp [-£ *(T) ■&}. 

Let us note here that the function, <£(£), contains com- 
plete information about the solution to the original prob- 
lem, Eq. (pQ), including the overall quantum phase accu- 
mulated by the wave-function during the time evolution 
(as we shall see below, this phase is of particular interest 
to the problem of dielectric response of TLSs in super- 
conducting qubits). An interesting question is whether 



and how this purely quantum phase can be restored from 
a solution of the corresponding classical Bloch equations 
that are usually considered in this context. Let us recall 
that a classical mapping can be achieved by introducing 
the average magnetic moment, 



m(i) 



* f (*)f*(*)- 



(6) 



Therefore, m 2 (t) = 1/4 and the classical equations 
of motion for the spin moment follow from dtm(t) = 



^(t) H(i),a #(£) and yield the familiar result 
d t m(t) = b(t) x m(t). 



(7) 



Let us recall that these Bloch equations are a saddle 
point of quantum spin dynamics, much in the same way 
that Newton's equations of motion governed by the force, 
[-VV(r)], represent a saddle point of the action describ- 
ing a quantum particle in the potential, V(r), and there- 
fore do not contain direct information about quantum 
interference and tunneling effects. Similarly, Eqs. ([7]) 
do not directly contain the quantum phase and to de- 
termine it one has to go back to the Schrodinger equa- 
tion. Another more abstract way to see this is by noticing 
that Eqs. ([7]) describe the motion on a two-dimensional 
(Bloch) sphere, m(t) G 5 2 , while the original quantum 
problem Eq. (|2j) describes motion on a three-dimensional 
sphere since S a h s (t) € SU(2) ~ S 3 . Now let us recall that 
there exists the Hopf fibration such that SU(2)/U(1) = 
5 2 , which summarizes the fact that classical equations, 
namely Eqs. ([7j), represent quantum motion modulo the 
U(l) phase dynamics. Fortunately, this phase dynamics 
can generally be restored from exact dependence of the 
m(t) solution, albeit in a non-local way. To see this, we 
can write the magnetization in terms of the ^-matrix as 

follows m(t) = ^ f (0) \&(t)&&(t)\ tf(0), where tf(0) 

and the corresponding m(0) = ^ (0)^^(0) are initial 
conditions for the wave-function and Bloch magnetiza- 
tion, correspondingly. Using again the well-known iden- 
tities for the Pauli matrices, we find the evolution matrix 
for the Bloch equations, m a (t) = R a p(i)mp(0), as fol- 
lows 

Rap(t) = Sap cos $ + n a np (1 — cos (j)) — e a p 7 n 7 sin <£>. 

(8) 
This three-dimensional matrix describes a rotation, 
R(t) G 50(3), and can be represented equivalently as 



R(t) = exp -&{t) • L 







where L = 




where L E so(3) ~ su(2) belong to the three-dimensional 
vector representation of the su(2) algebra. They are re- 

lated to the "usual" spin-1 representation (where J z is 
diagonal) via a simple linear transform. 

Therefore, we see that if we know an arbitrary solu- 
tion to the Bloch equations, m(£) we can at least in 



principle restore the function, <£(£), [see, Eqs. (J9j) and 
(JU], which uniquely determines the entire quantum so- 
lution. It also suggests that if we choose an arbitrary 
dynamic function on a sphere we may be able to restore 
the quantum Hamiltonian that would give rise to it, via 
mappings m(t) — » R(t) — » <&(t) -> S(t) — » H. However, 
the second step in this chain of transforms involves ef- 
fectively calculating a logarithm of the rotation matrix, 
which due to a complicated "analytic" structure of this 
matrix-logarithm function requires a careful calculation 
non-local in time. 

The subsequent Sections are devoted to constructing 
exactly solvable periodic-in-time Hamiltonians based on 
a specific Ansatz for the classical Bloch "magnetization," 
m(t). It further involves a restoration of the correspond- 
ing quantum U(l) phase via a straightforward integra- 
tion. More specifically, we "reverse engineer" the follow- 
ing Hamiltonian 



H = A t a x + f(t)a z 



(10) 



where f(t) = f(t + Tf) is a periodic function, with an 
arbitrary period, Tf. Our solution below also allows tun- 
ing of the average splitting, e = (f(t)) T _, and the AC 
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field amplitude, Af 



[fit) 



As mentioned 



in the introduction, this problem is of great importance 
to the physical problem of externally-driven TLS dynam- 
ics in superconducting qubits (with At corresponding to 
tunneling between the wells, e to a splitting of energy lev- 
els in a double- well potential, and Tf and Af being the 
period and the amplitude of the AC-electric field corre- 
spondingly). 

Our "guess" for the relevant Ansatz for the Bloch 
"magnetization," m(t), is based on a set of formal solu- 
tions discovered in the related problem of quenched dy- 
namics of fermionic superfluids.— ~ 21 i 24 i 25 Formally, the 
quenched dynamics of each individual Cooper pair is de- 
scribed by the Bogoliubov-de Gennes Hamiltonian, which 
is essentially a spin Hamiltonian that reduces to ([TO]) af- 
ter the unitary transformation a x — >• o z and a z — > —<J X , 
with At corresponding to a single particle energy level 
and f{t) to the superfluid order parameter. A realization 
of each particular form of the superfluid order parame- 
ter dynamics in a steady state can be unambiguously 
determined by the initial conditions 21 using the exact in- 
tegrability of BCS models Note that a self-consistency 
condition for the order parameter provides a limitation 
on the set of functions for which the corresponding prob- 
lem is integrable and for some initial conditions periodic- 
in-time self-consistent dynamics, /(£), can be realized. 
While in our TLS problem, there is no natural self- 
consistency constraint, such insights and constraints from 
the BCS problem help us narrow down the range of possi- 
ble Ansatze to restore reasonable physical Hamiltonians, 
which are also exactly solvable by construction. 

In what follows, we generalize the soliton analysis of 
Yuzbashyan 2 ^ and find a general soliton configuration, 
characterized by three independent parameters, which 



we denote as A± and A a . For the physical problem 
of interest, this conveniently implies that some, gener- 
ally speaking, non-trivial combination of these parame- 
ters will determine the arbitrary frequency, amplitude, 
and the dc-component of the field. Due to the periodic- 
ity, we can generally represent the AC-perturbation as a 
Fourier series 



f(t) = e + Af22f n cos(nuft). 



(11) 



(fT5|l can be cast to a more symmetric form, using another 
set of parameters A a and A±, which are chosen to be 
positive and are related to coefficients Cj as follows: 



c\ 



C2 



C3 



Aa 



(A 2 +-A 2 _), 






A 2 _+2A 2 ), 
A 2 )(A 2 -A 2 _ 



(16) 



Note that for certain specific choices of the parameters 
A±, a , the leading coefficient j\ ^> f n (n = 2,3,...) and 
one recovers the limit of a monochromatic AC-field, al- 
beit in the regime of weak driving (Affi <C max{A t , s}). 
Therefore, our non-linear analysis contains the standard 
linear response results as a simple special case. 



III. NON-DISSIPATIVE DYNAMICS OF THE 
AC-DRIVEN TLS 

In this Section we provide the details on the derivation 
of the exact solution for the TLS dynamics. We devote 
the special attention to the analysis of the U(l) phase of 
the wave function. We also elucidate the relations be- 
tween the parameters of our solution and the amplitude, 
phase and the dc-component of the external field, which 
may be useful for experimental applications of our theory. 



A. Ansatz 

We now focus on the Schrodinger equation for the half- 
integer spin in the magnetic field, h(t) = 2 (At, 0, f(t)). 
When written in terms of spinor components, it has the 
form 



f #+ = A t ^_ + f{t)1>+ 

\ iijj_ = A t V+ - JW- ' 

The corresponding Bloch equation is 

m(*) = 2(A t ,0,/(*))xm(t). 



(12) 



(13) 



Let us now make the following Ansatz for its exact solu- 
tion^ 

m x =D-Cf 2 , m y = Bf, m z = Af(t)+F. (14) 

From two of the Eqs. ^ we find A = 2A t £ and B = 
C. Thus among five parameters in ([T4]) only three are 
independent: F, B and D. The equation for the external 
field, /(£), can be obtained from ([T4]) using the condition 
m 2 = 1/4. This resulting equation for the function f(t) 
acquires the form 



/ 2 = -/ 4 -4c 2 / 2 + 8 Cl /-4c 3 , 



(15) 



where coefficients Cj are given by some combinations of 
parameters B,D and F [see Eqs. ([30]) below]. Equation 



Without loss of generality and to be more specific we 
also assume A + > A_ for the remainder of this paper, 
while A a can be assigned an arbitrary value. By virtue 
of expressions ([T6]) equation (fT5|) now reads 



/ 2 = [(/-A a ) 2 -A 2 _][A 2 h -(/ + A a ) 2 ]. 



(17) 



Below we will make several transformations that allow us 
to reduce (fT7|) to an equation for the Weierstrass elliptic 
function*^ Firstly, let us introduce a function, y(t), 



/(*) = A + 



2 
W) 



1 



A„ 



which satisfies the following equation 



4(y-a + )(y-a-)(y-l), x 



A+t 



(18) 



(19) 



where a± = 2A + /(A + + 2A a ± A_). Now, Eq. JE 
can be easily reduced to a well-known equation for the 
Weierstrass elliptic function by rescaling the parameters 
via the transformation 



so that 



y(x) = Z[x) 



§) -«*■ 



o+ + a_ + 1 



ei)(Z- e 2 )(Z-e 3 ), 



(20) 



(21) 



where parameters Cj satisfy the following conditions e\> 
62 > 63 and ei + e2 + e3 = 0. Coefficients Cj are de- 
termined by the parameters A a and A±. The specific 
expressions for the coefficients ej, however, depend on 
the relative values of the initially introduced set of pa- 
rameters and are given in Appendix A. Solution of the 
equation ([2T]) is 



Z(x) =V(x + uj'), J 






e 3 



(22) 



where V{x) is a Weierstrass elliptic function, K is 
a complete elliptic integral of the first kind 26 and 
k 1 = \/( e i — e 2)/(ei — es). Function Z{x) is a doubly- 
periodic function with the period along the physical time 
axis determined by, I — 2lj, where uo = K(^)/ v /ei — e% 
and k — y/l — k' 2 is a modulus of elliptic functions. Com- 
bining ([22]) with Eqs. ([20]) and ([18]) allows us to express 



fit) in terms of elliptic functions. Expression for fit) can 
be compactly written in terms of Jacobi elliptic functions. 
Just as it is the case for the parameters ej , the particular 
form of the resulting expression depends on the relation 
between A a and A± (see Appendix A). 



or a single soliton, where the period of the elliptic func- 
tion is taken to be infinite. Fig. 3 shows dynamic trajec- 
tories of the "magnetization" on the Bloch sphere given 
by exact Eq. ([T4]) that correspond to these particular 
/(^-dependencies. 
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FIG. 2: Plots of the function f(t) J23} in units of A+ for 
different values A a : (a) A a = 0.1 A+, A_ = 0.3A+; (b) 
A a = 0.5A+, A_ = 0.3A+; (c) A a = 0.3A+, A_ = 0.1A+; 
(d) A a = 0.5A+, A_ = 0.001 A+. We note that for the choice 
of the parameters (d) the period of f(t) diverges. The curves 
above are plotted for the value of At = 0.5A+. 

All cases considered here are summarized by the fol- 
lowing compact expression for the function, /(£), written 
in terms of Jacobi elliptic function sn as follows: 



where variable z is 
(t-t ) 



77 + sn 2 (z,ft) — 1 
77_sn 2 (z,ft) + 1 



A„, 



(23) 



[(A + +2A a )2-A 2 _]( ei -e 3 ) (24) 



and to = — cj / v /a + a_/A + . If we consider A± fixed, then 
the parameters rj± are given by one of the following ex- 
pressions depending on the value of A a : 



V± 



ei-e 3 

1 



±i, 

±K 2 , 



A a > 



A++A_ 



A+-A_ 



< A„, < 



A++A_ 



(25) 



ei— e 3 ' 



A a < 



Fig. 2 displays some representative dependencies of the 
driving field from the class of solutions described by 
Eqs. d23]), ([24]), and (|25j). Note that the curve in 
Fig. 2a is visually indistinguishable from a harmonic pe- 
riodic signal, Fig. lb and Fig. lc contain apparent non- 
monochromatic contributions to the periodic signal, and 
finally Fig. 2d provides an example of a degenerate case, 
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FIG. 3: TLS dynamics on the Bloch sphere (|13I14[) . Trajec- 
tories of TLS for the solutions described by Eq. ()23j) and 
depicted in Fig. 2 for the various set of parameters A a and 
A±. The latter take the same values used on Fig. 2. 



From the expression for the external field ([23]) it is, 
however, not immediately clear what set of parameters 
correspond to the regimes of weak and strong ac-driving. 
To clarify this issue, let us re- write ([23]) in the form more 
useful for practical applications. Let us first explicitly 
derive the amplitude, frequency and the dc-component 
of function fit). The period and the amplitude of oscil- 
lations of fit) can be immediately deduced from (j23|24p : 



4K(«) 



[(A + + 2A a )2_A 2 _](e 1 -e3) 



Af = ^ 



A + (r]++V- 



(26) 



V- 



1 



Lastly, the average value of the function fit) over its 
period is 



</(*)> 



A+?7 + 



1 _(^ ± +^) n( _ 7? _ )K) 



A = e, 
(27) 



with K(^) and IJ(n, k) being an complete elliptic inte- 
gral of the first and third kind correspondingly. As we 
have already mentioned, quantity ([27]) describes the dc- 
component of the external field. One can view Eqs. ([26] 
[27]) as the definition of yet another set of parameters A f , 
ujf = 27r/Xy and e = (/(t)), which allows us to cast ex- 
ternal field fit) into the form given by ([TT]) . We plot the 



dependence of these parameters on the ratio A_/A + in 
Fig. 3] for different values of A a while keeping the value 
of A t fixed. As we can see from [4] the limits of strong and 
weak ac-driving are easily attainable with the frame of 
our solution. In particular, we see that the regime of the 
strong ac-driving should be achieved for moderate values 
of A a and A_/A + -0.2. 
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FIG. 4: Plots of the amplitude Af, frequency ujf and dc- 
component e of the external field /(£), Eq. ([TT]) : (a) A a = 
0.1A+, At = 0.3A+; (b) A a = 0.5A+, A t = 0.3A+. 



Expressions ()23|24|25|) constitute one of the main tech- 
nical results of this paper. To get a further insight into 
the properties of our solution we refer the reader to Ap- 
pendix B where we consider few limiting cases for the 
function ([23]) . Quite generally, our solution represents 
the superposition of monochromatic waves with frequen- 
cies integer multiples of Uf = 2ir/Tf. As discussed in 
the Appendix B, solution ([23]) can be reduced to the 
monochromatic wave with frequency 2A + when A a = 
and A_ ~ A + . 



B. Wave function 

Having determined the form of the periodic field f(t) 
we employ the relations (J6j) to compute the amplitudes 
ip+(t) and ip-(t). First let us represent these functions 
as follows^ 



^±(t) = hM*)|e 



Ti<t>(t) Ja(t) 



(28) 



From these expressions, it follows that absolute values 
of the components ^ + and ip- as well as their relative 
phase (j){t) are determined by the instantaneous value of 
magnetization ([11]) . From Eqs. (|6|14p . we find 



i^±(*)i 

tan[20(t)] 



± 2A t Bf(t) ± F, 
f 



(29) 



(D/B)-p(tY 
where parameters B and D are determined from 



| = 2(A t 2 -c 2 ), B 
F = -2 Cl B/A t . 



1 



4J(A t 2 -c 2 p + 



A? 



c 3 (30) 



and parameters Cj's are given by ([16]) . Note that appar- 
ent ambiguity in signs for the parameters B and D as well 
as for parameter F is resolved by fulfilling the condition 



,2 _ 



1/4. 



C. Restoring the U(l) phase 



It has been mentioned above that the common phase 
a(t) has to be determined from the solution of the equa- 
tions ([T2]) . At first sight the resulting equation for a(t) 
appears to be very complicated, but it can be significantly 
simplified using Eqs. ([29]) . so that 



(31) 



. = 1 f(t)m z (t) 
° 2[l/4 -ml(t)Y 

After some algebraic manipulations, we find 



a(t) = I { A t 

o 
F 



4 



di 



[f 2 {f)-d\ p(t')-di 



(32) 



2S 2 [P{t') - d\] [p(t')-cP_ 



dt + ao, 



where a$ is determined by the initial conditions, d\ = 
(1/2B) ± D/B. One can evaluate the integral in ([32]) 
exactly and express in terms of elliptic a and ( functions 
(see Apendix C for details of this calculation). We note 
that on the grounds of Floquet theory we can represent 
an expression for the phase a(t) as a sum of two terms: 



a(t) - a = -i/£ + 7(£), 



(33) 



where j(t) = j(t + Tf) is a periodic function and v is a 
constant. Analytic expression for both of these quantities 
can be extracted from the analytic expression for a(t) 
listed in Appendix C. For example, from ([33]) it follows 
v = [a(t) - a(t + Tf)]/T f . In the limit when A a = 
and A + = A_ we find v = (A? + A^_) 1/2 , while in 
the limit when A a = and A_ = we obtain v = 



A t . For a general set of values A a and A± the resulting 
expression for v is not as simple as those listed above. 
For practical purposes, however, one can construct an 
approximate expression for v. By analyzing the behavior 
of a(t) numerically we find that for A a = 0, frequency v 
can be approximated (see Fig. [5^) by: 



Tf 



v(A a = 0) 



L f 



A t 2 



\t)dt. 



(34) 



We find qualitatively different behavior of v as a func- 
tion of 5 = A_/A + for nonzero A a . In that case, there 
appears to be a discontinuity in v at some critical value 
of A_/A + . The source of this discontinuity at least for 
small A a lies in the fact that a(t) oc m z (t) changes sign 
during its time evolution. For non-zero A a there are al- 
ways exists 5 C such that m z (t = 0) = 0, while for S > S c 
one observes m z (0 < t\^ < Tf) = 0. The sign change in 
m z (t) implies that the derivative of the quantum phase 
will change sign also (J3T]h so that the subsequent inte- 
gration yields the value of v smaller than the one found 
for 5 < 5 C , Fig. EJo. 
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FIG. 5: Dependence of exponent v as a function of the ratio 
A_/A+ for various values of A a . On panel (a) we compare 
the result of numerical computation of v from (|32|) and com- 
pare them with approximate expression ([34]) when A a = 0. 
Panel (b) shows the dependence of v for A a / 0. 

In order to get further insight into the physical mean- 
ing of the quantity v, we can employ the analogy be- 
tween the TLS and spin- 1/2 and define the magnetization 



M a (i) = (Vg(t)\a a \Wg(t))/2, where ^ g (t) is a general so- 
lution of the Schrodinger equation and can be expressed 
as a linear combination of the particular solution Sfr(t) 
(see below). Then one can show 2 ^ that the dynamics of 
the vector M can be represented as a linear superposi- 
tion of vector rh(t) precessing with the frequency of the 
field f(t) and a vector h(t) such that h • m = 0. Each 
component of the latter oscillates with frequency v. Our 
results from Fig. 3b suggest that the rate of precession 
of vector h will be significantly reduced as one tunes the 
parameter A_/A + . 

The solution of the Schrodinger equation we described 
above is only a particular solution from which the general 
solution can be constructed straightforwardly by taking 
advantage of the underlying symmetries of Eqs. ([T2]) . A 
general solution for the wave function, Sfr^ = (V^, ?A-)? 
can be presented as 



* fl (t) = Ci 






c 2 






(35) 



where Ci^ are integration constants, which satisfy 
|Ci| 2 + |C 2 1 2 = 1 and are to be determined from the 
initial conditions. For example, for the specific choice of 
an initial condition when the TLS at t = resides in one 
of its two states, 



* 9 (0) 



the coefficients Ci^ are 

Ci = v;(o), c 2 = v-(o). 



(36) 



(37) 



Expressions listed in this subsection amount to full de- 
scription of the ac-driven dynamics of an isolated TLS. 
In the next Section, we will briefly outline several appli- 
cations of our theory. For simplicity, we will mostly focus 
on the properties of the non-dissipative dynamics. 



IV. EXPERIMENTAL MANIFESTATIONS 

In this Section we discuss the physical behavior of sev- 
eral quantities which can be probed experimentally for 
various physical realizations of the TLS. Before we pro- 
ceed with the discussion on the application of our results 
and computation of physical observables, we derive the 
expression for the evolution operator and the density ma- 
trix which will allow us to compute probabilities which 
characterize the dynamics of the TLS. 

Evolution operator S(t) is defined by 



* fl (t) = 5(t)* fl (0). 



(38) 



From expressions ([35]) one can always write down a gen- 
eral expression for the evolution operator, which is valid 
for arbitrary initial conditions: 



S(t) 









iMo) 



-V+(o) 



(39) 
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Note that it is now straightforward to derive the density 
matrix from ([39]) using the following expression^ 



p(t) = S(t)p S\t), 



(40) 



where po is the density matrix of an initial state of the 
TLS. The expressions (|39|40p can be used as a basis to 
analyzed the effects of the environment dissipation on 
the dynamics of the TLS. In particular, one can deter- 
mine the probability of the TLS to remain in the initially 
prepared state P^^(t). 



A. Coherent destruction of tunneling 

The phenomenon of the coherent destruction of tunnel- 
ing (CDT) has been predicted theoretically^— for vari- 
ous physical realizations. Qualitatively, this phenomenon 
can be interpreted as the dynamical trapping of the TLS 
in one its states. For example, CDT occurs when the 
survival probability of the initial state dynamically ap- 
proaches unity. This phenomenon has its counterpart 
known in literature as driving induced tunneling oscilla- 
tions. This effect has been first analyzed theoretically in 
a series of papers^"— and observed experimentally for 
the first time by Nakamura et alM To compute the sur- 
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FIG. 6: Plots of the return probability Pt-+t(*)» E< 3 (SD> 
in the limit of the strong ac-driving: (a) e = 8.5At, Af = 
13.5A t , Uf = 21A t ; (b) e = 0.1A t , Af = 13.5A t , u f = 8A t . 



([40]) . It is, however, easier to use an expression for the 
wave function ([35]) with the initial conditions (|36|37p . In 
particular, let us choose the initial amplitudes such that 
both C\ and C2 are real and introduce an angle #, so 
that d = cos(i?/2). 

After some algebra, we obtain the following expression 
for the return probability 



Pt^tW =0 + cos ^ ' m *(*) + 



■ sin $ cos[2a(t)] • \j — — m 2 z (t) . 



(41) 



The CDT occurs when P^_^(£) « 1 and we assume the 
initial conditions ([36]) . From ([4T]) it follows that if we 
perform an averaging over time frame longer than Tf 
and T/j, = 27r/^, the third term in ([^T]) averages out to 
zero, so that employing ([13]) we find 



(P Ht (t))~- + 2cosM t B e 



ci 



(42) 



This equation approximately determines the parameter 
range for which CDT occurs. Fig. 5 displays representa- 
tive results for the return probability and illustrates the 
CDT phenomenon: As we can see, in the limit of strong 
driving, i.e. when Af ^> A t and cjf ^> A t , the return 
probability remains of order unity, which implies that 
the tunneling processes become strongly suppressed. We 
also have found that CDT remains robust and is present 
as long as the parameters A a and A± are such that 
the dynamics of the TLSs remains in the strong driv- 
ing regime. This qualitative behaviour of tunneling was 
found to be essentially independent of the ratio e/A t . 
These our findings agree qualitatively with the results 
reported previously in Ref. [15| for the monochromatic 
AC-field. Finally, we note that if a system of charged 
TLSs, e.g. OH-rotors present in AI2O3 dielectrics, is 
driven into such non-linear CDT regime by an external 
AC electric field, then the TLS tunneling and the cor- 
responding dipole polarization dynamics will be strongly 
reduced. This suggest that a strong non-linear drive may 
actually correspond to lower losses. 



B. Dielectric response 

Fig. 1 provides a pictorial example of a TLS charged 
defect, - an OH-rotor, which is one of the most likely 
candidates of physical two-level-systems responsible for 
dielectric losses in superconducting qubits. This rotor 
has a non-zero dipole moment p and, therefore, responds 
to an applied external electric field £(t). In the absence 
of other interactions which may affect TLS dynamics, the 
Hamiltonian describing the dynamics is ([TO]) with f(t) = 
e +p • £{t). By construction, the average dipole moment 
of an isolated TLS is determined by the following average 
within the spin mapping 5 



vival probability P^_^(£) we can use the density matrix 



d(t) = m z (t)p. 



(43) 
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The linear dielectric response function can be computed 
from (pl3]) by differentiating the corresponding compo- 
nents of the average dipole moment with respect to the 
amplitude of an external field £. To define a non-linear 
dielectric response function corresponding to a solution 
m z (t), which generally is a complicated function of the 
amplitude, we consider the spin-spin correlation function. 
Up to a pre- factor, given by the angle between the initial 
direction of the dipole moment relative to the external 
electric field, the dielectric response of an isolated TLS 
is defined by the Fourier component of the following cor- 
relator 5 



oo 

e{uj)= l -je i " t e- t 'ma z {t),d z {0)])dt 



(44) 



where square brackets denote a commutator between the 
corresponding spin operators. The exponential prefactor 
describes the dissipative effects of the environment and 
averaging is taken over the initial state of the TLS. We 
are introducing the dissipative effects on a phenomeno- 
logical level only and ignore the difference between the 
relaxation and dephasing processes. This is sufficient to 
get insight into the general properties of the exact spec- 
trum of the dielectric response due to an ensemble of 
identical TLS. Operators a a (t) in Eq. ([44]) correspond to 
the Heisenberg representation: 



& z (t) = S\t)aJ(t). 



(45) 



We remind the reader that formally the evolution opera- 
tor S(t) is given by 



S(t) =Texp 



z 

i f H(t')dt' 



(46) 



with T being a time-ordering operator. Using ([39]) and 
assuming the initial conditions \^(0) = (a*, 6*), for the 
correlator JC(t) = j{[&z(f), <J Z (0)]) under the integral in 
(|44]) we find: 



K{i) = 4(|a| 2 - \b\ 2 )m ± (t) sin[2a(t)] + 8ilm[a*b] ) 
x {m z (0)m x (t) — m x (0)m±(t) cos[2a(t)]}. 



(47) 



Here m±(t) = \A/4 — m l(t) an d we have fixed the ini- 
tial value of the field so that /(0) = 0. The subsequent 
time integration yields an expression for the dielectric re- 
sponse function. Analytic analysis of the response func- 
tion e(uj) is hindered by the fact that the correlation func- 
tion JC(t) is only a quasi-periodic function of time, since 
it is expressed as a combination of two periodic functions 
with different periods Tf and Th ([33]) . so that we have 
to resort to numerical calculation. In Fig. 6, we present 
representative plots of the real and imaginary part of 
e(uj) (|H|) for the initial conditions ([36]) . To interpret our 
results, we recall that the common phase, a(£), can be 
written as a sum of a linear-in-t term plus a periodic 
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FIG. 7: Plot of real and imaginary part of the response 
function e(cj). Note that discontinuities in real part and the 
peaks in imaginary part of the response function appear at 
frequencies cjdis = 2(ncj zb u), (n = 0, ±1, ±2...) in agreement 
with expression (|48|) . These plots has been obtained for the 
following values of the parameters: A a = 0.15A+, A_ = 
0.3 A+ and A t = 0.5A+. 



function ([33]) . Since m±(i) is a periodic function with 
the period Tj, we can express the corresponding terms 
in ([47]) in a Fourier series. Subsequent time integration 
yields a response function of the following type: 



«m= £ 



Xr, 



2nuj f ±2v+ Y 



(48) 



where e n are the corresponding Fourier coefficients. From 
this expression, we see that the peaks in the imaginary 
part of the response function describing the energy losses 
due to TLSs appear at frequencies, uj^is — 2(nuj ± v), 
commensurate with the driving frequency but with an 
overall shift determined by the quantum mechanical 
phase collected by a TLS over one cycle, v. Note that 
within the linear response theory, one typically keeps only 
the lowest Fourier harmonic in the spectrum (n = 0) and 
neglects all others. For the case of a monochromatic field 
the Fourier component with n = is kept so that 



Qin(^) = £0 



1 



^ '.»-L2ai/+ i -' 



a=± 



UJ ■ 



(49) 



and we recover the textbook result for the dielectric re- 
sponse function^ For a fixed set of parameters, however, 
one would only keep the largest contribution to the imag- 
inary part of e(uj). For example, the imaginary part of 
e(uu) is the largest for uu* c± 2v. Note also that apart from 
a difference in the value of the resonant frequency (which 
in the regime of weak driving is given by the energy that 
governs a stationary time evolution of a TLS eigenstate 
in the absence of any perturbations, v = yA? + e<2 i ^ ne 
response functions for the case of monochromatic field 
and the field given by ([23]) are qualitatively the same. 

In the array of non- interacting TLS, the response func- 
tion must be averaged over a distribution of the barrier 
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heights, direction of the electric field etc. We leave the 
detailed analysis directly applicable to the array of non- 
interacting and pairwise interacting TLS for a future pub- 
lication. 



V. CONCLUSIONS 

In this paper we presented an exact solution for the 
problem of AC-driven dynamics of a generic two-level 
system. Our approach was based on constructing a non- 
linear differential equation for the driving field, which has 
admitted an exact solution. The key feature of our so- 
lution for the external field is that it is fully described 
by three independent parameters. We have shown that 
one can interpret different nonlinear combinations of 
these parameters as an amplitude, frequency and a DC- 
component of the field. Being very general in nature, we 
believe that our results and methods can be applied to 
a wide variety of experiments ranging from NMR to the 
analysis of dielectric losses in amorphous materials. 

This research was supported by the Intelligence Ad- 
vanced Research Projects Activity (IARPA) through the 
US Army Research Office award W911NF-09-1-0351. 
Authors would like to thank Slava Dobrovitski and Ro- 
man Lutchyn for discussions related to this work. 



Appendix A: calculation of the parameters ej 

In this Appendix, we provide explicit expressions for 
the parameters e/s, which determine the explicit form 
of our exact solution for the external field ([22]) . As men- 
tioned in the main text, the particular expressions for 
these parameters, ej, depend on the relative values of 
A a and A±. For the choice corresponding to 



A a > 



A_+A + 



(Al) 



we have 



=(i) 



i 



(2 — a + — a_), 



1 



4 = §( 2a - ~ a + _1 )' 



(A2) 



=(i) 



1 



e 3~' = g( 2a +-a- - 1 )- 



In the opposite case of 

A+-A_ 



< A < 



A_ + A+ 



(A3) 



we have 



,(2) 



o( 2 ) 



-{2a- -a+ -1), 
-(2 — a_ — a+), 



(A4) 



4 2) = -(2a+-a_-l). 



Finally when 



we have 



A„< 



A+-A_ 



(A5) 



4 3) = -(2a_-a+-l), 
4 3) = l(2o+-o_-l), 



(A6) 



= (3) 



-(2 — a_ — a + ). 



We also remind the reader, that the coefficients a± in the 
above equations are given by a± = 2A + /(A + + 2A a ± 
A_). 



Appendix B: Exact solution for the function /(£): 
special cases 

In this Appendix, we consider a few special cases, 
where the exact solution given by ([23]) (which is gen- 
erally described by three independent parameters) is re- 
duced to a degenerate function with simpler properties, 
which is characterized by two parameters only. The first 
case we consider corresponds to A a — >• 0. As shown be- 
low, the choice of A_ c± A + , corresponds to an external 
field of the following form, f(t) ~ A+[l + qcos(2A+t)} 
with g« 1. Another case considered here is the limit 
A_ — >• 0, but with both A a and A + kept finite. In that 
case, f(t) can be represented as a single isolated soliton. 



1. Limit of A a -► 

Our goal here is to recover the limiting case for our 
solution corresponding to S a = 0. It can be shown that 
in this limit, 



l-k' 



k' 



y/l - k' 2 



1 + lfe'' 

and the solution for the driving field reads: 

2 



/(*) 



|p[fcA + 



(t-to)] + l-e 3 



■}■ 



(Bl) 



(B2) 



We demonstrate below that the expression in the brackets 
can be cast into single Jacobi elliptic function dn(A + £, k). 
For this, we use the identity 



P ' / 

such that it gives 

/(*)=A+ 



e 3 + (ei - e 3 ) 



1 



sn 2 (i£,ft) ' 



sn 2 (ii,^) - (ei - e 3 ) 



sn 2 (^, k) + (ei - e 3 )_ 



(B3) 



(B4) 
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and variable u equals u = ^(1 + 5-)A+t. This expres- 
sion can be further simplified by means of the following 
relation between the Jacobi elliptic functions: 



dn(iii, fti) 



1 — ftsn 2 ^, k) 
1 + ftsn 2 (u, k) 



where 



u\ — (1 + ft)w, ACi = 



2y^ 
l + « 



Indeed, from expressions (|A6|) for A a = we have 
e 2 - e 3 1 - (J_ 1 



(B5) 



(B6) 



(B7) 



ei - e 3 1 + 5- ei - e 3 
so that Ki = k and we find: 

/(*) = -A + dn[A+(*-t ),fc]. (B8) 

Finally, when fc -)► (A_ -> A+) it follows^ that 

f(t) ~ -A+ [1 + gcos(2A + t)] , g < 1, (B9) 

We find that for the special values of parameters the line 
shape of the external field is given by the cosine. 



According to the expressions above for that case k = 
and sn(>,0) = sin(,/4A2 - A\t/2). It follows: 



/ = -A a - 



4A + (2A a -A + )/(2A a + A + ) 



l-cos(^4A2-A2.t)+2 



2A a -A + 
2A a +A + 



(B13) 



We see that when A_ = external field has a line shape 
of a single pulse. Note that our solutions (JB11)B13|) do 
not contradict to our assumption of the periodicity of f(t) 
since both these solutions correspond to the case where 
the period of f(t) goes to infinity. 

Appendix C: calculation of the common phase a(t) 



In this Appendix, we outline the main steps, which 
allow to compute the integral ([32]) exactly. The cal- 
culation includes the following transform of the special 
functions involved that reduces the integrand to a form 
amenable for exact integration of the Weierstrass elliptic 
function: 26 



2. Limit A_ 







To derive an explicit form of the driving field, /(£), in 
this case, we work with the general solution (|23j) . Let us 
first assume that 

A a < A+/2. 

T Then, the case A_ = corresponds to k — 1, which in 
turn implies 



sn(/u, 1) = tanh(u), u = -Xt 



(BIO) 



and A = -v/A^_ — 4A^. After some simple algebra, we 
find 



Aa' 



A 2 



2A a - A+ cosh(At) 



(Bll) 



which up to the minus sign, is exactly the same expression 
as the listed in Ref. [25]. Finally, let us consider the 
parameter range with 



A a > A+/2. 



(B12) 



/ 



aV(u) 



<yV(u) 



6 a 

— du =—u 



log 



aS — /3j 
<j(u + v) 



a(u — v) 



2u((v) 



(CI) 



where a,(3,j,5 are some constants, a parameter v is de- 
termined from the derivative of the Weierstass function, 
V{v) = —(5/7, cr(u), and ((u) are the Weierstrass elliptic 
sigma and zeta functions.— 

The next step is to write down the function, /(£), ex- 
plicitly in terms of the Weierstrass function. Combining 
expressions ([18120122)) , we have: 



/(*) = -A+ 



V(X + LJ')-1 



V(x + io') + l-ej 



A a 



A + t 

yoqro 



(C2) 



and 



where a± = 2A + /(A + + 2A a ± A_), x 

j = 1, 2 or 3 depending on the value of A a (see Appendix 
A). Let us now consider the first integral in 
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I 



d\ 



-dt 



j _ V a + a - 



d+^^ + wO + l-ej] 



P{t')-d\^ v 2A+ 7 \(A + + A a -d + )P(x / + a;0 + (Aa-d+)(l-e i )-A + (l + e i ) ( ^ + 



r 



-d_|_) > eta/ 
(C3) 



Here the index j of the coefficient ej is determined by 
the value of A a (see Appendix A). The remaining terms 
can be written in a similar form and the corresponding 
integrals can be evaluated using (|C1|) . as we have done 



for the first one (|C3|) . Since the resulting expressions for 
the a(t) turn out to be too cumbersome, we do not list 
them here. 
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